function EEG = cvc_load_bdf(bdffile,fs,elpfile)
% CVC_LOAD_BDF
%
%  Synopsis
%  ========
%
%  EEG = cvc_load_bdf(bdffile,fs)
%  EEG = cvc_load_bdf(bdffile,fs,chanlocsfile)
%
%  -- Author: Mads Dyrholm --
%     Center for Visual Cognition, University of Copenhagen.
%     2009 - March 2010
%
%  Purpose
%  =======
%
%  Read Biosemi BDF file and create EEGLAB EEG-structure 
%  with sample rate decimation. Zero-phase anti-aliasing
%  lowpass filtering is used.
%
%  Inputs
%  ======
%
%  bdffile - Filename of BDF file.
%              
%  fs - New sample rate. Must be integer fraction of 
%  BDF sample rate, or it will be rounded to nearest 
%  integer fraction (at CVC the BDF sample rate is 2048Hz).
% 
%  chanlocsfile - File that comes out of using the 
%  'cvc256.seq' sequence with the digitizer software
%  (.elp or .sfp). Use 'cvc_standard256.elp' for lookup 
%  locations with clamped fiducials.
%
%  Outputs
%  =======
%
%  EEG - EEGLAB structure. It has a custom '.cvc' field
%  which will be used by CVC_PREPARE_SPM.

if (nargin<3) elpfile = []; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       maxrecs = 50;        % if you get OUT OF MEMORY, reduce this number !
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

edf = openbdf(bdffile);
trigchan = edf.Head.NS
breakflag = (edf.Head.NRec > maxrecs);
if breakflag==1,
  BLOCKRECORDS = maxrecs;
else
  BLOCKRECORDS = edf.Head.NRec;
end
ALLEEG = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% decimation filter design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
edfFs = edf.Head.SampleRate(1);
r = edfFs/fs;
if r~=round(r)
  r = round(r); % decimation
  fs = edFs/r;  
  fprintf('\n\n Warning!!!: Output sample rate adjusted to %f\n\n',fs);
end
nfilt = 8;
rip = 0.05;
[b,a] = cheby1(nfilt, rip, .8/r);
u = zeros(edfFs/2+1,1); u(round(length(u)/2))=1;
h=filtfilt(b,a,u);
delay = round(length(u)/2);
memory = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% getting the channel locations
if ~isempty(elpfile)
  fprintf('Getting the channel locations.\n');
  EEG_chans = pop_chanedit([],  'load', elpfile ); %, 'filetype', 'polhemus'});
  %EEG.cvc.elpchanlocs = EEG_chans;
  % extra fiducial dummys
  verylast = length(EEG_chans);
  %try
  %  EEG_fids = EEG_chans.chaninfo.nodatchans;
  %  for fidx=1:3
  %    EEG_fids(fidx).X = EEG_chans.chanlocs(verylast-3+fidx).X;   
  %    EEG_fids(fidx).Y = EEG_chans.chanlocs(verylast-3+fidx).Y;
  %    EEG_fids(fidx).Z = EEG_chans.chanlocs(verylast-3+fidx).Z;
  %    EEG_fids(fidx).sph_radius = sqrt(EEG_fids(fidx).X^2 + EEG_fids(fidx).Y^2 + EEG_fids(fidx).Z^2); % we dont use this, but anyway
  %  end 
  %  EEG.cvc.elpfidlocs = EEG_fids;
  %end
  
  EEG_chans = pop_chanedit(EEG_chans, 'delete', 1); % ref
  for dummy = 1:3 % remove weird fiducial channels
    EEG_chans = pop_chanedit(EEG_chans, 'delete', 1);
  end
  for dummy = length(EEG_chans)+(0:-1:-2) % FID extra fiducial dummy channels
    EEG_chans = pop_chanedit(EEG_chans, 'delete', dummy);
    %EEG_chans = pop_chanedit(EEG_chans, 'changefield', { dummy, 'type', 'FID'}, 'changefield', {dummy, 'datachan', 0 });
  end
else 
  elfile = []; 
end

EEG_chans

% get events, and event channel
fprintf('Getting the events.\n');
EEG_events = pop_biosig16cat(bdffile, 'channels', trigchan, 'rmeventchan', 'off');%, 'channels', datachannels, 'ref', refchannels);
%EEG_events = pop_biosig(bdffile, 'channels', trigchan, 'rmeventchan', 'off');%, 'channels', datachannels, 'ref', refchannels);

[types,numbers] = eeg_eventtypes(EEG_events)

% allocate
downsampleddata = zeros(edf.Head.NS-1,ceil(EEG_events.pnts/r),'single');

% running indexes
dnidx = 1;
downphase = 0;

% get busy
for recordblock = 0:BLOCKRECORDS:edf.Head.NRec
  fprintf('\n********************\nProgress: %i%%\n********************\n',round(100*recordblock/edf.Head.NRec));
  
  % read block of data
  RECORDS = [recordblock min(edf.Head.NRec+1,recordblock+BLOCKRECORDS)]; 
  EEG = pop_biosig(bdffile, 'blockrange', RECORDS); %, 'rmeventchan', 'off', 'channels', datachannels), 'ref', refchannels);
  EEG.data = double(EEG.data);
  
  % for mysterious reason pop_biosig might give us the transposed data
  [M,N] = size(EEG.data);
  if (M>N), EEG.data = EEG.data'; end %... so fix it!
  pnts = max(M,N);
  
  % correct somewhat for huge DC offsets
  if (recordblock==0) 
    chanoffset = mean(EEG.data,2);
  end
  for ch=1:size(EEG.data,1)
    EEG.data(ch,:) = EEG.data(ch,:) - chanoffset(ch);
  end
  
  % filter.
  fprintf('Anti-aliasing filter w. constant group delay and downsampling.\n');
  [EEG.data,memory] = filter(h,1,EEG.data,memory,2);
  
  % downsample
  tmp = downsample(EEG.data',r,downphase)';
  
  % correct for filter delay
  if (recordblock==0)
    tmp = tmp(:,max(round(delay/r),1):end);
  end
  
  % insert
  updidx = size(tmp,2);
  downsampleddata(:,dnidx:dnidx+updidx-1) = single(tmp);
  dnidx = dnidx + updidx;
  
  % prepare for next downsampling
  downphase = mod(pnts-downphase,r);
  if (downphase>0) downphase = r-downphase; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

fprintf('\n Decimating the event latencies . . !\n');

for ev = 1:length(EEG_events.event) 
  EEG_events.event(ev).latency = max((EEG_events.event(ev).latency-1) / r +1, 1);
  EEG_events.urevent(ev).latency = max((EEG_events.urevent(ev).latency-1) / r +1, 1);
end

EEG_ref = EEG.ref;
EEG = eeg_emptyset;
EEG.data = downsampleddata;
EEG.srate = fs;
EEG.ref = EEG_ref;
EEG.event = EEG_events.event;
EEG.urevent = EEG_events.urevent;
EEG = eeg_checkset(EEG);
if ~isempty(elpfile),
  nEEGdigitized = length(EEG_chans);
  % empty dummy channels
  
  %  % EEGLAB 8
  %  comextra = {};
  %  for nn=nEEGdigitized+1:size(EEG.data,1) 
  %    comextra = cat(2, comextra, {'append', nn-1, 'changefield' , { nn ,  'labels', sprintf('e%i',nn)} ,  'changefield' , { nn , 'type', 'EEG'},  'changefield' , { nn ,  'datachan', 1 } });
  %  end
  EEG.chanlocs = EEG_chans;
  %  
  %  keyboard
  %  EEG=pop_chanedit(EEG, comextra{:});

  for nn=1:nEEGdigitized
    EEG.chanlocs(nn).type = 'EEG';
  end
  
  for nn=nEEGdigitized+1:size(EEG.data,1) 
    EEG.chanlocs(end+1).labels = ['e' num2str(nn)];
  end
  
  EEG = eeg_checkset(EEG);
  % fit locations to sphere
  EEG.chanlocs = pop_chancenter(EEG.chanlocs, [], nEEGdigitized+1:size(EEG.data,1));
  EEG = eeg_checkset(EEG);
end
EEG = eeg_checkset(EEG);
fprintf('\n DONE. \n');
